//
// Created by Ivor_ on 2023/3/4.
//

#include "ordinaryDifferentialEquation.h"

double* eulerExplicitOde(int numOfPoints, double* variableVals, double functionValueInit, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValueInit;
    for (int i = 1; i < numOfPoints; ++i) {
        functionValues[i] = functionValues[i-1] + derivativeFunction(variableVals[i-1], functionValues[i-1])*(variableVals[i] - variableVals[i-1]);
    }
    return functionValues;
}

double* eulerImplicitIterationOde(int numOfPoints, double* variableVals, double functionValueInit, double accuracy, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValueInit;
    for (int i = 1; i < numOfPoints; ++i) {
        double iterVariable = functionValues[i-1] + derivativeFunction(variableVals[i-1], functionValues[i-1])*(variableVals[i] - variableVals[i-1]);
        double iterVariableTemp;
        do {
            iterVariableTemp = iterVariable;
            iterVariable = functionValues[i-1] + derivativeFunction(variableVals[i], iterVariable)*(variableVals[i] - variableVals[i-1]);
        } while (fabs(iterVariable - iterVariableTemp) > accuracy);
        functionValues[i] = functionValues[i-1] + derivativeFunction(variableVals[i], iterVariable)*(variableVals[i] - variableVals[i-1]);
    }
    return functionValues;
}

double* eulerImplicitForecastCorrectOde(int numOfPoints, double* variableVals, double functionValueInit, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValueInit;
    for (int i = 1; i < numOfPoints; ++i) {
        double forecast = functionValues[i-1] + derivativeFunction(variableVals[i-1], functionValues[i-1])*(variableVals[i] - variableVals[i-1]);
        functionValues[i] = functionValues[i-1] + derivativeFunction(variableVals[i], forecast)*(variableVals[i] - variableVals[i-1]);
    }
    return functionValues;
}

double* rungeKutta2Ode(int numOfPoints, double* variableVals, double functionValueInit, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValueInit;
    for (int i = 1; i < numOfPoints; ++i) {
        double k1 = derivativeFunction(variableVals[i-1], functionValues[i-1]);
        double k2 = derivativeFunction(variableVals[i], functionValues[i-1] + (variableVals[i] - variableVals[i-1])*k1);
        functionValues[i] = functionValues[i-1] + .5*(k1 + k2)*(variableVals[i] - variableVals[i-1]);
    }
    return functionValues;
}

double* rungeKutta3Ode(int numOfPoints, double* variableVals, double functionValueInit, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValueInit;
    for (int i = 1; i < numOfPoints; ++i) {
        double k1 = derivativeFunction(variableVals[i-1], functionValues[i-1]);
        double k2 = derivativeFunction(.5*(variableVals[i-1] + variableVals[i]), functionValues[i-1] + .5*(variableVals[i] - variableVals[i-1])*k1);
        double k3 = derivativeFunction(variableVals[i], functionValues[i-1] + (variableVals[i] - variableVals[i-1])*(-k1 + 2.*k2));
        functionValues[i] = functionValues[i-1] + (k1 + 4.*k2 + k3)*(variableVals[i] - variableVals[i-1])/6.;
    }
    return functionValues;
}

double* rungeKutta4Ode(int numOfPoints, double* variableVals, double functionValueInit, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValueInit;
    for (int i = 1; i < numOfPoints; ++i) {
        double k1 = derivativeFunction(variableVals[i-1], functionValues[i-1]);
        double k2 = derivativeFunction(.5*(variableVals[i-1] + variableVals[i]), functionValues[i-1] + .5*(variableVals[i] - variableVals[i-1])*k1);
        double k3 = derivativeFunction(.5*(variableVals[i-1] + variableVals[i]), functionValues[i-1] + .5*(variableVals[i] - variableVals[i-1])*k2);
        double k4 = derivativeFunction(variableVals[i], functionValues[i-1] + (variableVals[i] - variableVals[i-1])*k3);
        functionValues[i] = functionValues[i-1] + (k1 + 2.*k2 + 2.*k3 + k4)*(variableVals[i] - variableVals[i-1])/6.;
    }
    return functionValues;
}

double* adamsExplict2Ode(int numOfPoints, double* variableVals, double* functionValuesInit, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValuesInit[0];
    functionValues[1] = functionValuesInit[1];
    double stepSize = variableVals[1] - variableVals[0];
    for (int i = 2; i < numOfPoints; ++i) {
        functionValues[i] = functionValues[i-1] + .5*stepSize*(3.*derivativeFunction(variableVals[i-2], functionValues[i-2]) - derivativeFunction(variableVals[i-1], functionValues[i-1]));
    }
    return functionValues;
}

double* adamsExplict3Ode(int numOfPoints, double* variableVals, double* functionValuesInit, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValuesInit[0];
    functionValues[1] = functionValuesInit[1];
    functionValues[2] = functionValuesInit[2];
    double stepSize = variableVals[1] - variableVals[0];
    for (int i = 3; i < numOfPoints; ++i) {
        functionValues[i] = functionValues[i-1] + stepSize/12.*(23.*derivativeFunction(variableVals[i-3], functionValues[i-3])
                - 16.*derivativeFunction(variableVals[i-2], functionValues[i-2]) + 5.*derivativeFunction(variableVals[i-1], functionValues[i-1]));
    }
    return functionValues;
}

double* adamsExplict4Ode(int numOfPoints, double* variableVals, double* functionValuesInit, double(*derivativeFunction)()) {
    double* functionValues = (double*)calloc(numOfPoints, sizeof(double));
    functionValues[0] = functionValuesInit[0];
    functionValues[1] = functionValuesInit[1];
    functionValues[2] = functionValuesInit[2];
    functionValues[3] = functionValuesInit[3];
    double stepSize = variableVals[1] - variableVals[0];
    for (int i = 4; i < numOfPoints; ++i) {
        functionValues[i] = functionValues[i-1] + stepSize/24.*(55.*derivativeFunction(variableVals[i-3], functionValues[i-3])
                - 59.*derivativeFunction(variableVals[i-3], functionValues[i-3]) + 37.*derivativeFunction(variableVals[i-2], functionValues[i-2])
                - 9.*derivativeFunction(variableVals[i-1], functionValues[i-1]));
    }
    return functionValues;
}
